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vM Abstract A micro-hydromechanical model for granular materials is presented. It 

5-H combines the discrete element method (DEM) for the modeling of the solid phase 

and a pore-scale finite volume (PFV) formulation for the fiow of an incompressible 
pore fiuid. The coupling equations are derived and contrasted against the equa- 
r~^ tions of conventional poroelasticity. An analogy is found between the DEM-PFV 

y—{ coupling and Blot's theory in the limit case of incompressible phases. The simu- 

lation of an oedometer test validates the coupling scheme and demonstrates the 
',^_,' ability of the model to capture strong poromechanical effects. A detailed analy- 

^^ sis of microscale strain and stress confirms the analogy with poroelasticity. An 

£Jq immersed deposition problem is finally simulated and shows the potential of the 

method to handle phase transitions. 



Keywords Discrete Element Method, Pore-Scale Finite Volumes, Poromechanics, 
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1 Introduction 



The description of the mechanics of saturated porous media is a problem of great 

interest in engineering. Various approaches have been attempted in order to model 

the behaviour of such systems, the complexity of their structure and the interaction 

^^ between the solid and the fluid (liquid, gas) phases. At the microscopic (sub-pore) 
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scale, the solid and fluid phases occupy different portions of the spatial domain 
and interact at their common interface. Thus, the microscopic fields which describe 
the properties of constituents may be considered as continua within a single phase, 
while exhibiting discontinuities at the interfaces between phases. A key feature of 
saturated porous media is the two-way coupling between the deformation of the 
solid matrix and the fluid pressure. This coupling is the source of the so-called 
poromechanical effects, which govern for instance the consolidation process. Karl 
von Terzaghi was the first author to describe these effects, in the framework of 
continuum mechanics [35]. Later, Maurice Biot adopted Karl von Terzaghi's ideas 
and extended the one-dimensional consolidation theory to the three-dimensional 
case for linear elastic porous solids [3]. 

The main objective of the work presented in this paper is to capture the 
poromechanical effects in a discrete numerical model of a granular material. The 
model will be formulated for the special case of incompressible phases. In this sit- 
uation Terzaghi's equations are consistent with Blot's theory of poroelasticity, and 
the poromechanical effects are maximized J^ . We will restrict the study to quasi- 
static behaviour. The choice of following a discrete approach in modeling granular 
materials became more and more popular in the last decades. Analyzing the be- 
haviour of granular materials at the particles scale enables a better insight into 
phenomena which are governed by particles behaviour, and gives access to kine- 
matic and static variables which are extremely difficult to measure in experiments. 
Hence the links between the microscale and the macroscale can be investigated. 
A popular method for discrete modeling is the discrete element method (DEM). 
Initially, the DEM has been developed without considering the effect of fluids, 
and was therefore restricted to dry granular materials. In the recent years, how- 
ever, many efforts have been devoted to couple the DEM with models of fluid flow 
within the medium. The various methods that have been developed differ in the 
modeling techniques adopted for the description of the fluid flow and are briefly 
reviewed below. 

Continuum-based models use continuum formulations and coarse-grid meshing 
for the fluid part of the problem |25ll33l[38lll2] . The coarse mesh defines subdo- 
mains of the DEM model over which the porosity and velocity of the solid phase 
are averaged and introduced as field variables of the continuum formulation. In 
turns, the solution of the continuum problem gives fluid velocity and momentum 
exchange between the phases at each node of the grid. The latter is discretized 
again with appropriate rules to determine forces to be applied on the particles 
of the DEM model. The solution can be obtained with classical numerical meth- 
ods such as finite differences (FD) [25] or finite volumes (FV) 03] • Continuum 
approaches generally lead to affordable CPU costs, since the number of degrees 
of freedom (DOFs) associated to the fiuid can be much smaller than the number 
of solid particles. They need a series of phenomenological assumptions and rely 
on empirical relations such as Ergun's relation [17 . A drawback is that they may 
require a calibration procedure for each new type of microstructure. More impor- 
tantly, they are inherently unable to describe accurately the effects of fiuid at the 
particles scale. This is because the variables of the fluid problem are averaged at 
the scale of the coarse-grid cells, whose size is generally larger than the average 
size of solid particles. 

Inversely, microscale models are based on a very fine discretization of the void 
space. A Navier-Stokes problem where the boundary conditions are prescribed at 
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the surface of each particle is then solved using conventional techniques (e.g. the 
finite element method (FEM) fTS^), or particle based methods such as the Lattice- 
Boltzmann method (LBM) 5,28,22,26]. The advantage of such direct methods 
is that, unlike continuum-based models, they do not rely on phenomenological 
assumptions. Typically, the fluid-grain interactions are governed solely by a no- 
slip condition at the fluid-solid interfaces. The limitations of these models come 
from the high computational cost, especially when applied to the solution of a 
nonlinear system of equations in three dimensions. The high costs are due to the 
need for a very small mesh size compared to the size of particles. Even if the 
LBM is known to minimize this computational cost, the minimal number of nodes 
necessary for accurate results in three dimensions still restricts the problem size to 
small numbers of particles [21]. In addition, it should be noted that incompressible 
flow represents a difficulty for the explicit methods such as conventional LBM- 
BGK, which are inherently based on density fluctuations. This feature makes the 
use of the explicit methods difhcult in the situations of strong poromechanical 
coupling associated to incompressible or very weakly compressible phases. 

To overcome the high computational cost of the microscale models, without 
introducing all the phenomenological assumptions of continuum-based methods, 
pore-network (PN) models are a very attractive compromise [24 . PN models are 
based on a representation of the void space as a network of connected pores and 
throats, where the properties of the throats are supposed to reflect the effect of lo- 
cal void geometry on the flow. This modeling technique greatly reduces the number 
of unknowns as compared to micro-scale models. PN models have been extensively 
used for studying single-phase and multiphase flow in rigid porous materials, but 
they have been rarely coupled to models of deformable solid skeleton, although 
some examples exist in two dimensions for granular materials [20111] and in three 
dimensions for fractured rocks [23] . A numerical model of the PN type has been 
developed recently for incompressible flow in sphere packings [10]. In this model, 
the spatial discretization leads to fluid elements whose sizes are of the same order 
as the size of the solid particles, thus introducing an intermediate scale between 
the microscale and the continuum scales mentioned above. At the same time it 
can reflect the effects of the local geometry of the pore space on the flow and 
gives predictive estimates of the permeability 36J , unlike continuum based meth- 
ods. Hereafter, we refer to this method as the Pore-scale Finite Volumes method 
(PFV). In the present paper, we present a coupling between the PFV method and 
a DEM model in three dimensions. 

In the first part, the governing equations of the DEM and PFV models are 
recalled, and the coupling equations are established and discussed. The analogy of 
the system of equations with Blot's equations is highlighted in the particular case 
of incompressible phases. By such analogy, the DEM-PFV coupled model should 
be able to recover results of classical poroelasticity in boundary values problems, 
provided they share similar assumptions. This is the case for the simulation of an 
oedometer test, which is presented in the third part. In the fourth and last part, 
the coupling is applied to the simulation of a granular deposition problem. We will 
introduce microscale definitions of the stress and strain tensors for detailed anal- 
ysis of the results. The evolution of stress, strain, and fluid pressure is presented 
and discussed. The onset of liquefaction events and fluid-solid transitions which 
are observed in the material behaviour, show the interest of following a discrete 
approach in simulating fluid-particle systems. 
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Fig. 1 Spherical particles interacting at contact points. 



2 Numerical model 



2.1 Solid phase 



2.1.1 Explicit DEM 



The discrete element method (DEM) has been extensively used to study soil and 
rock mechanics, providing, for instance, some insights into shear strength and 
deformation properties of geomaterials. The approach is fully micromechanical, the 
solid phase being modelled by defining the mechanical properties of the interaction 
between the grains that compose it 13 . In what follows, we recall the generic 
aspects of the DEM method for the simulation of systems of spheres. The algorithm 
presented uses the explicit finite difference scheme for time integration and assumes 
smooth contact behaviour. Those two aspects are sufficient requirements of our 
fiuid-coupling algorithm. A large majority of publicly available softwares (including 
commercial and open-source) are using similar schemes according to our survey. It 
offers a good fiexibility for incorporating additional physical effects and couplings 
(see e.g. [32]). Hence, the coupling scheme we are developing can be potentially 
combined with a large number of existing DEM codes. Implementation details 
of the open-source code Yade [3l], used for the present study, can be found in 
Smilauer and Chareyre [37] . 

We consider a system of A'' spheres interacting with each other at contacts (fig. 
[I]). The kinematics of each sphere is described by six degrees of freedom (DOFs). 
We note X^ = {xi,0i} the generalized position of sphere i, with Xi the position 
of its center of mass and 0, the rotation represented as a R'^ vector (although in 
Yade rotations are represented by quaternions, a rotation vector is mathematically 
equivalent and more convenient in the present derivation). The symbol X used 
without indice will refer to a vector containing all the DOFs of the system (6 x N) 
and will be called the global position vector. Similarly, x will contain all the 
translational DOFs. 
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The translational motion of each sphere in the system is governed by Newton's 
second law, which relates the forces exerted on a particle to its acceleration: 



miXi = / a''nds+ / p^^dv, (1) 

where rrii is the mass of a particle occupying the volume 71 , cr^ n is the stress 
applied at the particle surface in the direction of the unit normal n, p^ is the mass 
density and g is the gravitational acceleration. The first integral is the total force 
exerted on particle i by the other particles, we note this force Ff . 

Assuming that the contact areas are negligible, the contact interactions can be 
represented by contact forces f/j^ acting at contact points: 

miXi = ^ fik + mi g. (2) 

fe=o 

where ric is the total number of contact points. 

By introducing the global vectors x = {xj} and F'^ = {Ff } containing respec- 
tively the positions and forces for the 3 x A'' translational DOFs in the system, the 
evolution can be defined by the matrix relation: 

x = M"^(F'= + W), (3) 

where M is the global mass matrix and the components of W contains the gravi- 
tational terms. 

A relation similar to eq|2] holds for the rotational DOFs by replacing force and 
mass by the torque and inertia tensor, respectively. Hence, the system of eq[3]must 
be supplemented with 3 x A' equations for the rotational DOFs for completeness, 
in order to define completely the dynamics of the system: 

X = J-^(T^-hW'), (4) 

with J the generalized inertia matrix, T'^ the generalized force vector (including 
torques) (W' is simply W with zeros appended for the rotational DOFs). 

The contact forces f/j^ appearing in eqJ2] are computed according to contact 
laws, which often describe irreversible behaviour. For this reason, it is not possible 
in general to define a unique relation between positions and forces. Instead, it is 
usually only possible to define the rate of change of f^ as a function of positions 
and their time derivatives: 

file = 6ifc(Xi,Xk,Xi,Xk) (5) 

where the function bi^ defines the constitutive behaviour of the contact between 
particles i and k. Then, in terms of the global components, we can introduce an 
operator B (nonlinear in general) corresponding to the summation of all forces 
and torques on the particles, so that 

F'^ = S(X,X) (6) 

The explicit DEM method consists in integrating equations |4] and [6] with a 
time-stepping algorithm, updating positions and forces at each step. 
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The most common algorithm is based on a discretization of acceleration (left 
hand side of eqlsl) with a centered second order finite difference scheme, which 
reads 

It results in an explicit equation for computing 'X.t+At- 

The contact forces are updated in the new configuration according to eqJ5J 
where Fj , ^^ /g and Xt-^-At/2 Suie replaced by second order approximations at mid- 
step: 

Ff+At = Ft + B{'X.t + At,Xt + At/2)'^t 

^t+At - Xt (8) 



Xt+zit/2 ^^ 

We may remark that this matrix representation of the problem is not very 
common in the literature on DEM. It is introduced here for an easier presentation 
and discussion of the coupled problem in the next sections. 

2.1.2 Elastic-frictional contacts 

In the simulations that follow, the contacts will be modelized using an elastic- 
plastic contact law 13J. It does not imply any loss of genericity of the DEM- 
PFV coupling algorithm, which only depends on the explicit integration method 
formerly described. 

The contact force is decomposed into its normal and tangential parts. The 
normal force (compressive only) is proportional to the normal displacement dn 
and to the normal stiffness kn, 

f _ j -knd„ if dn < , . 

■^" " \ if d„ > ^^> 

dn is defined as d„ = \\xi — Xj\\ — Rt — Rj, where Ri and Rj are the radii of 
the particles in contact, so that dn = when the spheres are exactly tangent. The 
shear force depends on the shear stiffness kt, and is integrated using eqjs] where the 
relative velocity dt at contact point Xc depends on the spin ujij of each particle. 

^•^"^*'^* . ^ ■ . ^ (10) 

dt = Xj -I- UJj X (Xc - Xj) - Xi - tJi X (Xc - Xi) 

Coulomb's friction is introduced by imposing a maximum magnitude for ft: 

||ft|| <tan(,?!.)/„ (11) 

where (f> is the angle of Coulomb's friction at contacts. 

kn and kt are defined, for each couple of particles i-j in contact, as functions 
of an elastic modulus E and a non-dimensional constant a = kt/kn- 

'^^ - ^ (fli+flj) (12) 

rZt ^ CI Kn 

This definition of stiffness is convenient for defining problems independently of the 
mean particles size, since it results in a constant ratio between E and the effective 
bulk modulus for a given type of packing. 
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2.2 Fluid phase 

2. 2. 1 Viscous flow 

We assume that the porous medium is saturated by an incompressible fluid whose 
flow is governed at the micro-scale by Stokes equations, which express fluid mass 
and moment conservation at small Reynolds and large Stokes numbers. Situations 
in which the Reynolds number is small are called viscous flows, because viscous 
forces arising from shearing motions of the fluid predominate over inertial forces 
associated with acceleration or deceleration of fluid particles. The Reynolds num- 
ber is defined as: 

p-^u-Vu 

Re = -^ 13 

where p-* is the fluid density, u its velocity, /i its viscosity. By noting as [7 a 
characteristic fluid velocity for the problem and d a characteristic length describing 
the problem (e.g. the average throat dimension), the Reynolds number can be 
computed as, 

Re=^''^^ (14) 

In fluid-particle systems, the Stokes number is a dimensionless parameter that is 
defined as the ratio between a viscous diffusion term and a term related to the 
fluid acceleration: 

St = ^ (15) 

P 'at 

Using the same notation of eq]14[ and noting as Tc a characteristic time of rear- 
rangement of the particles, Stokes number can be written as: 

The motion of the particles produces in fact a poral flow whose relaxation time is 
of the order of magnitude of r^, = (P /i^ (remember that i/ = fi/p-^). A situation in 
which Ti, « Tc means that the flow induced by the particles displacement attains 
rapidly a new equilibrium. In this case {St » 1), the acceleration of the fluid can 
be neglected and the hypothesis of steady laminar flow holds. 

For small Reynolds and large Stokes numbers, Stokes equations, for the de- 
scription of the flow, read: 

Vp = pV^u (17) 

V-u = (18) 

where u, and p are the microscopic fluid velocity and piezometric pressure, re- 
spectively. The piezometric pressure p is related to the absolute pressure p" via 
p = p"" — p* gz, with g the gravitational acceleration and z the depth coordinate. 
A no-slip condition is assumed for the fluid velocity at the grain boundaries. 

Doing the hypothesis of slow viscous flow, the flow and the forces exerted on 
the particles can be computed using the PFV model recently developed [10] . The 
key aspects of this method are summarized below. 
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Fig. 2 Regular triangulation (A) and its dual Voronoi diagram (B). 



2.2. Z The PFV model 

Regular triangulation and its dual Voronoi graph are used to discretize the void 
space. A system of tetrahedra arises from the triangulation in a 3D framework, 
each tethraedron representing a pore (see figJ2|A) and figl3|. The vertices of the 
triangulations are the spheres of the packing, so that the displacement of the 
particles will be reflected by the deformations of the mesh elements. The dual 
Voronoi diagram constitutes a network whose edges do not cross the solid phase 
(see figpFB)). Such network represents the flow path within the porous sample and 
allows the formulation and resolution of the flow problem, practically upscaling 
Stokes equations at the pores scale. The pressure field is defined by the values 
of pressure in each tetrahedral element, located at the vertices of the Voronoi's 
graph. The mass balance equation integrated on the pore i, and recast into a 
surface integral by using the divergence theorem, gives: 



V/ 



/ (u-v 



) • nds = ^ Qi 



(19) 



where V/ is the total pore volume, (u — v) is the velocity of the fluid relative to 
that of the solid phase. The integral on facet Sij gives a fiux exchanged between 
adjacent tetrahedra, noted qij (see figisp). The sum of fiuxes equals the rate of 

volume change of the pore V/ , which is related to the velocity of the particles. 
This sum can be seen as a discrete divergence operator applied on fluid velocity. 



A key aspect of the model is the expression of the flux qij through a facet as 
a function of the local geometry and pressures in the pores. Stokes equations im- 
ply a linear relation between pressure gradients and fluxes. The expression can 
therefore take the form 

Hi 



Qij — Qij - 



(20) 



where a local pressure gradient is deflned as the ratio between the pressure drop 
Pi — pj and the distance lij between pores i and j (euclidean distance between the 
Voronoi vertices associated to each pore), gij is a term expressing the hydraulic 



conductance of the domain 0ij represented on figjSp. Combining equations 19 
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Fig. 3 Volume of fluid in a pore (A), adjacent pores and local connections (B), fluid domain 
of pore contour (C), pore partition for hydraulic radius definition (D) .10. 



and [20] gives a relation linking the discrete pressure field to the velocities of the 
particles: 



V, 



J=3l 



The following definition has been proposed and validated for gij [101136] : 



5»j 



of r>h ' 

2^ 



(21) 



(22) 



where Rij is the hydraulic radius (defined hereafter), Sf, is the area occupied by 
the fiuid in facet Sij , and /z is the fluid viscosity. The hydraulic radius was defined 
as the ratio between the volume &ij occupied by the fluid and the area 7^^ of 
solid-fluid interface (see flglsp). 



2.2.3 Forces 

The total force F; exerted by the fluid on a particle i results from the pressure 
and viscous stress acting at the surface: 



F[ = / (-p"n + Tn)ds 



(23) 



where dFi denotes the solid surface of the particle i, p°' the absolute pressure and 
T the viscous shear stress tensor. Remembering that p = p° — P gz, F; can be 
evaluated by summing three terms C|10jl: 



; = / —p^gznds+ / pnds + / t nds 
JaTi -Jdri -Jari 



(24) 



where the flrst term denotes the buoyancy force, the second term the integral of 
piezometric pressure, and the third term the integral of the viscous stress. The 
second and third term vanish if the fluid is at hydrostatic equilibrium (constant 
P)- 
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2.3 Hydromechanical coupling 

The previous studies on the PFV method were limited to flow in rigid assemblies 
of spheres. Therefore, the hydromechanical coupling was not considered. In this 
section, we detail the micro-scale coupling equations that appear when the DEM 
and PFV models are combined in a unified framework and we discuss the relation 
with the field equations of classical poroelasticity. 

2.3.1 Coupling equations 

The coupling is defined by two matricial relations. The first one corresponds to 
the mass conservation equation [2l] written for all pores, complemented with the 
boundary conditions: 

GP = Ex + Q, + Qp, (25) 

where G is the conductivity matrix containing terms gij/hj of equation |21[ P 
the column vector containing all values of pressure, and E is the matrix defining 

the rates of volume change such that V/ = (Ex)i. Q, and Qp are flux vectors 
reflecting boundary conditions, respectively source terms (imposed fluxes) in Qg 
and imposed pressures in Qp. Solving this equation gives the field of fiuid pressure 
P as function of particles velocity x. 

The second relation results from the addition of the fiuid forces to the contact 
forces in Newton's equation [3J which becomes: 

Mx = F'' + W + F-^, (26) 

where F'^ and F-* are the contributions of contact forces (as they appear in eqlsk 



and fluid forces (eq 24), respectively 



The expression of F-* of eq 24 can be expressed in matricial form as a function 
of the pressure field P and a matrix S refiecting the local geometry of the sphere 
packing (please refer to iO_ for details): 

F-^ = SP (27) 

2.3.2 Coupled problem 

Combining equations |25| and |26[ we end up with an explicit ordinary differential 
equation of order 2, where X is the only remaining unknown: 

X = M-'(F" + W + SG-'(E • X + Q, + Qp)). (28) 

Integrating this equation numerically poses no major theoretical difficulty. The 
methods used in finite element solvers may also be used here (see e.g. [15]). How- 
ever, attention must be paid to the method choosen in order to minimize the 
computational cost. For this purpose, we developed a semi-implicit scheme [71|8] 
which makes the resolution of the coupled problem possible with quite accept- 
able computation times. The overhead of running a coupled simulation with this 
scheme is of the order of 100% of the computation time of the same simulation 
without fiuid. 
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We note that the constitutive relations wiU to lead to energy dissipation in the 
simulated systems by either frictional effects at contacts (eq. 11) or viscous effects 
due to the fluid (eq. 22). It is relatively common in DEM to introduce additional 
sources of dissipation, damping the equation of motion (eq. [2]) artifially [11] for a 
faster convergence to static equilibrium. Such numerical damping is not required in 
the DEM-PFV since it is naturally damped, and it would in fact lead to erroneous 
results regarding the time evolution of the systems. In what follows, the results 
are obtained without introducing any form of numerical damping. 



2.4 Relation with classical poromechanics 

Although the coupling equations were obtained only from micro-scale considera- 
tions, we will show that they may be seen as a discrete form of the field equations 
of conventional Blot's theory of poroelasticity for quasi-static deformations. It is 
worth noting this feature since the next section presents the comparisons with 
Terzaghi's analytical solution of the monodimensional consolidation problem. 
In the case of incompressible phases (Biots coefficient = 1), the coupling equa- 
tions of Blot's theory are the Poisson's equation that result from the continuity of 
Darcy's velocity [15] : 

^(V-dx) + V-(-i^dVp)=0, (29) 

and the equation of local equilibrium: 

-V-(T' + Vp={l-n){p' -pf)g, (30) 

where p = p°' — p* g z \a the excess ("piezometric") pore pressure and Kd is the 
hydraulic conductivity of the medium, a' represents Terzaghi's effective stress and 
dx is the displacement of the solid phase. 

In small strain linear elasticity, cr' = C | ( Vdx -|- Vdx ) with C the stiffness 
tensor. Substituting cr' by this expression in eq |30| gives a Navier-type equation 
where Vp can be seen as a body force: 

-V ■ C^(Vdx + Vdx"^) -h Vp = (1 - n){p'' - /)g. (31) 

In the light of this system of partial differential equations, we can reconsider the 
equations of the pore-scale DEM-PFV formulation. Fi rstly , we observe that in the 



case of steady rates of deformation, X vanishes and eq 28 becomes an equilibrium 
equation. Secondly, we note that in the special case when all contacts behave 
purely elastically, a linear relation exists between the contact force vector F"^ and 
the generalized displacement dX via a global stiffness matrix C [1]. Hence, the 
system satisfies a relation of the form: 

C^ dX -h W -F SP = (32) 

Considering only one particle i of the system, interacting with ric particles in 



contact, and with nj incident fluid cells, eq 32 implies 



Y, Cfk(dX, - dXfc) + y. p^ g + Y. ^-kPk = 0. (33) 



fc=0 fc=0 
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In this equation Cffc is the so called rigidity matrix [1], that is multiplied by the 
relative displacements dXi — dXfc to give the contact forces; the second term 
expresses the weigth of particle i {Vi = particle volume); the last term sums the 
contributions of fluid pressure in incident pores. Equation [33] can be seen as the 
discrete form of a Navier equation, where X^fcZg' C^fe(- — dXfc) is an operator 
defined for the discrete displacement field, and equivalent to the operator Jy V • 
C{V • +V •) that would arise in the conventional finite volume formulation for 
a continuum. Similarly, the third term is a discrete operator equivalent to the 
integral of the pressure gradient /y V-. We note that the microscale counterpart 
of the divergence of the effective stress is a sum of contact forces. 

The analogy between eq |21| and eq ]29| is direct since in the former, by definition, 
the two terms represent the local rate of volume change of the pore space and the 
divergence of the fluid velocity averaged in the pore. 

Consistently, we also remark that the assembled matrices for the coupled 
boundary value problem, as given by equations [25] and ]32] do not differ from 
the ones obtained through the discretization of Blot's equations using the FVM 
[29| or FEM [15] methods. As a last note, we can remark that the sets of elements 
that we obtain from the triangulation and the tessellation systems are also quite 
similar to the ones found in unstructured FV mesh [271130] . 

From this comparison, we can conclude that the DEM-PFV coupled model may 
recover results of classical poroelasticity in boundary values problems, provided 
they share similar assumptions (this is further discussed in the next section). The 
ID diffusion problem known as Terzaghi's consolidation is well suited for such 
comparison. It is used as a benchmark test in the next section for the validation 
of the model. 



3 Oedometer test 

The consolidation process is a classical hydro-mechanical problem, of primary 
importance in geomechanics. In the case of incompressible phases and when the 
deformation occurs in only one direction, as in oedometer tests (fig. ]4|, Blot's 
theory of poroelasticity coincides with the classical solution of Terzaghi. 

The equation of monodimensional consolidation is a diffusion equation on p, 
that reads: 

^ =6* — 
dt dz^ 



C.^2 (34) 



where z is the coordinate along the axis of application of the load (vertical direc- 
tion), Cv the consolidation coefficient, defined as follows: 

Kd Eoed 



G^, 



p^ g 



(35) 



where K^ is the coefficient of hydraulic conductivity of the soil, Eoed = Aa^j As^ 
the oedometric modulus. A non-dimensional time parameter is introduced, Ti,, 
defined as: 



T^ = % (36) 



where t is the time, and Ld the longest drainage path for the generic fiuid particle. 
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By noting Z = z/ H, with H the height of the sample, the analytical solution 



of the problem of eq 34 (Taylor, 1948) , expressing the evolution of the interstitial 
pressure p, reads: 

P(^'^'') = E ,f\n • -"(| • (2™ + 1) • ^) • e-(f •(^-+^))^-" (37) 
^-^ n ■ (2m + 1) 2 

m=0 ^ ' 

As the interstitial pressure dissipates, the deformation of the sample takes 
place. The average degree of consolidation expresses the evolution of the total set- 
tlement over the final settlement that would result at the end of the consolidation 
process, and is defined as; 

^(^) -1 V ^ -{^■(2m + l))-.T, ,„gs 

Sfrrral Z^ (yr • (2m + 1))^ ^ ' 



3.1 Numerical simulations 

The ID consolidation problem will be used as a benchmark test for the DEM- 
PFV coupling. The numerical solution will be compared directly to the analytical 
solution of eqpT] and eq|38[ The DEM-PFV model is consistent with Terzaghi's 
assumptions under the following conditions: 

— cl) The deformation is small and the stress-strain behaviour is linear; 

— c2) Blot's coefficient = 1 (incompressible phases); 

— c3) The displacement is unidirectional; 

— c4) The permeability and the elastic properties are constant in space and time. 

The cl) condition can be obtained with the DEM for small deformations, as 
shown below. The c2) hypothesis holds in the DEM-PFV model since both phases 
are incompressible. Boundary conditions can be specified to be consistent with 
the c3) condition, but the displacement inside the material will not be purely 
unidirectional due to local fiuctuations in the displacement field. By definition, 
a discrete model generally does not satisfy the c4) condition, since the heteroge- 
neous arrangement of particles leads to fabric heterogeneities, then to a spatial 
variability of the elastic properties and the local permeabilities. In addition, the 
evolution of the microstructure during the deformation may induce an evolution 
of the mechanical and hydraulic properties of the medium. 

In a first simulation, however, the permeability and the oedometric modulus will 
be nearly constant in time since the deformation will be very small. In a second 
simulation, conditions cl) and c4) conditions will be relaxed and the consequences 
will be commented. 

The oedometer test was simulated on a cubic sample, sized 0.1m, bounded by 
rigid plates. The boundary conditions are shown on figure[4] Lateral displacements 
are prevented {uxx = Uyy = 0). The fluid pressure p = OkPa is imposed at 2; = 



and z = H (corresponding to a two-way drainage with Ld = H/2 in eq 36 1 . 
5000 slightly polydispersed grains were employed to build the sample. They were 
first compacted isotropically using the REFD growth algorithm ,9j, by which a 
confinement stress of 5 kPa is applied. A relatively dense sample was created 
(n ~ 0.36), to minimize the dispersion of pores' dimension and avoid strong spatial 
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Fig. 4 Ocdometric boundary conditions and DEM model of oedometer test (with less particles 
than in the actual simulation). 



Pressure 
5000. 




Fig. 5 Pressure field at T„ = 0.1. 



heterogeneities within the sample. Then, the oedometer test itself was simulated 
by applying an increment of stress Aoext = 5 kPa on the top plate. Table IT] gives 
the main parameters of the simulation. 

In addition to the oedometer test, two independent simulations were performed 
on the same sample. First, a compression test on the same sample but in dry 
conditions, with the same loading path as the oedometer test, was used to evaluate 
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Fig. 6 Simulated load-unload cycle on the dry sample with Aa^xt = 5 kPa (a) and AcText 
200 kPa (b). 



Table 1 Oedometer test - Input data of the test of fig.s 
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INPUT DATA 



Number of grains 


H 


5000 


pf 


[kg/m3] 


1000 


P" 


[kg/m^] 


2600 


Sample dimensions 


H 


0.66 X 0.66 X 0.66 


M 


[kPa ■ s] 


0.25 


c'so 


[m] 


0.0395 


Confinement stress 


[kPa] 


5 


Auext 


[kPa] 


5 


PO 


[kPa] 





E 


[kPa] 


15000 


kt/kn 


H 


0.5 



the oedometric modulus of the sample. In a loading-unloading cycle at Aa^xt = 
5 kPa, the response is completely reversible, showing that we are in the purely 
elastic regime at such stress-strain level (see fig|6]ja)). The final strain obtained 



with the dry sample (8.525-f 0~ ) defines an oedometric modulus Eoed = 5895 kPa. 
Second, a permeability test was simulated on the initial geometry by imposing p = 
at z = and p = 1 at z = H. The solution of this simulation gave fluxes that were 
used to determine the equivalent permeability of the sample. The permeability 
obtained from the permeameter simulation is Kd = 7.07623 • 10~^ m/s. 

The numerical results that were obtained for the oedometer test are summa- 
rized in table [2] FigJT] shows the evolution of excess pore pressure at half the 
total height of the sample. The plotted value is an average computed on the plane 
{x,y,z = H/2). Figis] (left) shows the evolution of pore pressure in space (z/H) 
and time (Tv). The pore pressure is normalized by the increment of applied stress 
Aaext- FigIs] (right) shows the evolution of the settlement s, normalized to the 
value s final obtained on the dry sample. 

The excess pore pressure (plotted on figJTl) rose up almost instantaneously to 
5 kPa (= aext = Pmax) and then gradually decreased. We remind that in Terzaghi's 
formulation, the initial condition is p = aext, at t = {Tv = 0). In our simulation. 
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Table 2 Ocdoineter test - Numerical result (see fig.s 
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RESULTS 



Final strain 

Eoed 

Kd 



6300 ■ 



5000 K 



ct 400C ' I 

■ — ■ I 1 



pi I 

i" 3000 J 



2300 4 



[kPa] 
[m/s] 



8.525 ■ 10^'' 

5895 

7.07623 ■ 10-5 



Time [s] 

Fig. 7 Oedometer test - Fluid pressure measured at 2 = H/2 during the consolidation process. 



p = initially, and the maximum value was obtained in finite time: p = a^xt 
at Tv = 3.2 • 10~^ ~ 0. This short time lag can be explained by the inertia 
of the system: pore pressure results from the velocity of the solid phase, which 
is not established instantaneously. This short delay is also found in experiments 
[19| and in other numerical simulations [5lll2j. For the rest of the process, the 
evolution of pressure in space and time is found to be in good agreement with the 
analytical solution, as it can be seen on figJ8](left). The same conclusion holds for 
the evolution of the settlement, as shown on figjs] (right). 



3.2 Microscopic stress and strain 

For a more detailed analysis of the results, the computation of microscale strain 
and stress tensors has been implemented in Yade-DEM code [34j . The tensors are 
defined in particle-centered volumes Ve and Va (fig. |9|. The microscopic stress 
associated to one particle is defined as a sum over the contacts [16l l2]: 



1 V^ c,fc 



Va 



(39) 
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Fig. 8 Oedoineter test - Evolution of pore pressure {Aa^xt = 5 kPa) (left) and settlement 
(right) 




Fig. 9 Particle-centered domains for the definition of micro-strain (a) and micro-stress (b). 



where x'^' is a contact point and f^' ' the corresponding force. Va is the reference 
volume associated to the particle in the Voronoi tesselation. Note that this tensor 
does not reflect the average stress in the solid material. For this purpose, we would 
have to divide by the volume of the particle, not by Vcy , and the stress applied by 
the fluid on the contour should be accounted for. Instead, this micro-stress only 
reflects that part of the external loading that acts through the contact network. 
For this reason, it can be seen as the microscale analogue of Terzaghi's effective 
stress, as will be confirmed by the results. 

The microscale strain tensor for one particle is defined as a function of the 
displacements of the particles adjacent to that particle in the regular triangulation, 
which define the polyhedral domain Ve- The average displacement gradient in an 
equivalent continuum that would be contained in Ve is 



< Vdx >= — / Vdxdw = 77 / 



dx ® nds. 



(40) 



where the displacement dX on the contour dVe is defined as a piecewise linear 
function, equal to the displacement of the particles at the vertices and linear on 
each facet. This expression generalizes in three dimension the expression proposed 
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Fig. 10 Consolidation problem - Tv = 0.10. On the left, profile of fluid pressure p, micro-stress 
a^z and micro strain Wzz- On the right, the same data are plotted on a a^z vs. Ezz graph. 
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in [14] . The micro-strain is then obtained as the symmetric part of the gradient; 



1 T 

-(< Vdx > -h < Vdx >-*) 



(41) 



Both tensors may not be very meaningful at the scale of one particle alone, 
and there is no clear constitutive relation between e and W at such a small scale 
[2]. They have however the interesting property to converge, respectively, to the 
average strain and stress in an equivalent continuum when they are averaged in 
larger domains containing many particles |6j . Hence, we can define stress and strain 
at a certain height z in the sample as weighted averages of the tensors associated 
to the particles present at (or near) this height. It let us plot profiles of stress and 
strain as functions of z. 

(left) shows the profiles of fluid pressure p, micro-strain Szz and 
at Tv ~ 0.10. As it was expected, fluid pressure and effective stress 



So, fig 



10 



micro-stress ffzz • 

profiles are complementary, as stated by the Terzaghi's effective stress principle 
a = a' + p. The linear relation between the micro-stress and the micro-strain can 

shows a SD-visualization 



be observed in the right diagram (fig 10 (right)). Fig 



11 



of the deformation field within the sample, at the same time Tv ~ 0.10. A profile 
similar to the ones of figJ5] can be observed, although the strain field reflects more 
local fluctuations than the pressure field. 



3.3 Nonlinear consolidation problem 

In this section, the influence of conditions cl) and c4) (see previous section) on 
the flnal solution is examined by simulating larger deformations. The input data 
are the same than the ones summarized in table [T] except that the amplitude of 
the stress increment is larger: Aaext = 200kPa. The flnal strain obtained through 
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Fig. 11 Strain field at Tv = 0.10 - 3D-visuaHzation. 



Table 3 Nonlinear oedometric consolidation - Numerical result. 



RESULTS 


Final strain 


[-] 


3.94-10-2 


Initial Permeability 


[m/s] 


7.07623e - 05 


Secant Permeability 


[m/s] 


6.41317e-05 


Initial -Boed 


[kPa] 


5895 


Secant E^ed 


[kPa] 


5620 


Initial Cv 


[-] 


0.0417 


Secant Cv 


[-] 


0.0360 



the dry compression test was 3.94 • 10~^. On figMright) it can be seen how the 
initial state was not recovered after unloading the sample. Table [3] summarizes the 
results for this new simulation. The initial and secant values of the permeability, 
and the intitial and secant values of the oedometric modulus are reported, as well 
as the respective consolidation coefficients, computed using eq |35[ 

In Fig|12[ the analytical solution is computed by using the initial (tangent) 
value of Cv Such choice leads to a discrepancy between the analytical and the 
numerical solutions. More precisely, the rate of deformation is underestimated. 
On the contrary, in fig]13[ the analytical solution is obtained with the secant 
value of Cv The agreement is rather good in this case, although the difference 
remains larger than in the linear case (fig. Isl). This is because in case of large 
deformations the permeability is no longer constant, and the stress-strain behavior 
is no longer linear. Well then, permeability changes and non-linearities discard the 
use of unique constants for the permeability and the oedometric modulus of the 
granular medium. The results obtained at large strain are thus picking out the 
limitation of the linear solution from Terzaghi. 
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Fig. 12 Ocdometric consolidation with large deformations - the analytical solution is com- 
puted with the initial (small-strain) value Cv: evolution of pore pressure (left) and settlement 

(right). 
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Fig. 13 Oedometric consolidation witli large deformations - the analytical solution is com- 
puted with the final (secant) value of Cv: evolution of pore pressure (left) and settlement 
(right). 



4 Immersed granular deposition 



In this section the evolution of the fluid pressure and the effective stress is observed 
in a granular deposition problem. The simulation was set in order to simulate a 
fluid-filled vessel in which a number of spheres is placed and left to deposit under 
the action of gravity. No fluxes are allowed through lateral and bottom boundaries. 
A cloud of immersed spheres is created, as represented on fig|I4| Table [4] reports 
the parameters of the simulation. The fluid pressure is recorded at six different 
heights during the simulation, according to the scheme represented on fig]I4| 
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Fig. 14 Simulation of an immersed granular deposition (a). Position of fluid pressure sensors 

V^ (b). 



Table 4 Immersed granular deposition - Input data. 



INPUT DATA 


Number of grains 


[-] 




5000 


Sample dimensions 


H 


2.00 


X 1.50 X 0.75 


M 


[kPa • s] 




0.10 


d^o 


[m] 




0.06 


P" 


[kg/m»] 




2600 


pf 


[kg/m3] 




1000 


Pext 


[kPa] 







E 


[kPa] 




15000 


kt/K 


H 




0.5 



4.1 Critical gradient 

Once the sedimentation process reaches a steady rate, the weight of the soil is 
completely carried by the fluid phase, that gets over-pressurized. These conditions 
are mechanically equivalent to a fluidization of the soil, which by definition consists 
in the cancellation of the effective stress as an effect of an upward ground water 
seepage. This situation occurs when the hydraulic gradient equalizes the so-called 
critical gradient 



Vp 

7/ 



7sat - 7/ 

7/ 



(42) 



where 7sat indicates the specific weight of the mixture: "fsat = 7s(l — n) + ^f{n), 
with 7s and 7^ the specific weight of the solid and the fluid phase, respectively. 
Considering the initial porosity of the suspension, n = 0.61, we can directly eval- 
uate the critical gradient of the simulated suspension: ic = 0.624. 
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Fig. 15 Immersed granular deposition, case of tabUI Fluid pressure measurements (a). Hy- 
draulic gradient (b). 



4.2 Simulation results 



The result in terms of fluid pressure measured along the height of the vessel is 



plotted on fig 15 a). The gradient of fluid pressure, evaluated for each couple of 
consecutive sensors, is plotted as well (b). On fig 15 [a), when the pressure measured 
at a certain layer equals the one measured at an adjacent layer, it means that at 
that depth the spheres are newly in contact and the layer has stabilized. It can 
be observed, on fig|15[b), how the initial values of the hydraulic gradient are close 
to the critical gradient defined previously, and finally goes to zero as the packing 
stabilizes. On fig]16|a number of stress states, characterizing the evolution of the 
simulation, is represented. The initial load is entirely carried by the fluid phase, 
and the effective stress is none (A, see also flg 15 for the meaning of the capital 



letters) . As soon as the particles start touching each other in the lower strata, the 
consolidation process starts, and the stress is transferred from the liquid phase to 
the solid skeleton (B-C). Finally, once the fluid pressure is fully dissipated and the 
consolidation process completed, the load is entirely carried by the solid skeleton 
(D). 

These results highlight the flexibility of the coupling, which can handle the 
poroelastic couplings but also the transitions between a loose suspension and a 
solid state. 



4-2.1 Stokes and Reynolds numbers 



As introduced in section |2.2[ the soundness of the steady laminar flow hypothe- 
sis, in case of fluid-particles systems, is evaluated through the evaluation of the 
Stokes number. It is important to give an estimation of the Stokes number for the 
simulation of the granular deposition, especially for the initial conflguration which 
is characterized by a relatively high value of porosity. From eq |16[ we may use 
the average distance between the particles (average throat diameter) as a charac- 
teristic length of the problem (d in the equation). Noting as e this distance and 



Title Suppressed Due to Excessive Length 



23 




■■ ■" Effective Stress 



EffettJve Stress | 
Fluid pressure 



K 



D 



\ 



lMi5 5000 5000 4000 500C OOCIO 7005 5000 



V Effective stress 
— Fluid pre^sLi^e 



ICOC 2000 50OO 4000 MOO OOOO 7000 3000 



1000 2000 5000 400C 50OO fOOO 7000 30OC 



Fig. 16 Effective stress and fluid pressure evolution during tfie granular deposit. 



knowing the initial value of porosity n = 0.61, we have: 



n= 1 



Vt 



4/3 7rdio 
8(d5o + e)3 



0.61 



e = 0.014m 



(43) 



where dso is the mass-median-diameter of particles, Vs and Vt are the solid and 
the total volume, respectively. Thus, we obtain: 



St 



500 >> 1 



(44) 



where the order of magnitude of Tc is estimated by looking at fig |15| how much 
time the lowest layers take to get stabilized, thus Tc ~ 30s. The Reynolds number 
can be also estimated by considering the permeability of the sample at the initial 
conditions, which is Kd = 4.01621 • 10~^ m/s, and the critical hydraulic gradient. 
We obtain Re ~ 10~^. The hypothesis of slow viscous flow holds for this simulation, 
thus validating the flow model assumptions. 



5 Conclusions 



The paper is devoted to the presentation of the DEM-PFV coupled model for 
viscous flow in granular materials. A key feature of the model is the description 
of the interaction between the solid and the fluid phases at the scale of pores 
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and particles. This strategy, inspired by the pore-network approach, significantly 
reduces the computational cost generally associated to the direct simulation of 
flow in porous materials. Although the model is formulated in the framework of 
discrete mechanics, it results in a system of equations that is consistent with the 
theory of poroelasticity. 

The ability of the model to solve transient problems in the quasi-static regime 
has been evaluated in the oedometer test simulation. The solution obtained, in the 
case of small deformations, is quantitatively in good agreement with Terzaghi's 
analytical solution, in terms of evolution of the excess pore pressure, stress and 
settlements, in time and space. In the case of large deformations, the limitations 
of Terzaghi's formulation due to the use of constant values of permeability and 
oedometric modulus have been highlighted. Micro-scale definitions of strain and 
stress have been introduced and allow meaningful interpretations of the results. 
Namely, the micro-scale stress can be considered as a micromechanical analogue 
of Terzaghi's effective stress (this may be also true to some extent in unsaturated 
materials, as suggested by [32]). 

To our knowledge, it is the first time that such quantitative agreement with the 
analytical solution of Terzaghi's problem is reported for a DEM-based hydrome- 
chanical model (although relatively good agreement was also found in \TTJ. The 
reason may be that the implicit PFV formulation is based on strict incompressibil- 
ity of the fiuid, while the other coupling methods often use explicit formulations 
which need a finite (strictly positive) value of compressibility. The time-step of an 
explicit method depends on fiuid compressibilty, so that capturing the porome- 
chanical effects involved in the incompressible limit requires extremely small time- 
steps. Since, in addition, the accuracy of conventional methods (CFD,LBM,SPH) 
requires thousands of fiuid DOFs per solid particle (while the number of PFV 
DOFs is of the order of the number of particles), the poromechanical effects may 
be recovered only at the price of inconveniently high computational costs. 

Since the 1-D consolidation problem of Terzaghi is one of the simplest boundary 
value problems governed by poromechanical effects, we believe that it could be a 
standard problem for benchmarking coupling models. It could be argued that 
Terzaghi's problem is assuming incompressibility for the fiuid while water, for 
instance, is not incompressible. However, the compressibility of water is so low 
that many experiments on soils are not infiuenced by it. Therefore, the numerical 
models should be able to produce results corresponding to the incompressible limit. 

The numerical result that was obtained by simulating the granular deposition 
problem confirmed the robustness of the model and its ability to handle a wide 
range of solid-fiuid interactions. Hence the strength of the DEM-PFV coupling is 
twofold. First, it is consistent with poromechanics although the granular material 
is defined only by its microscale geometry and contact laws. Second, the model 
is not restricted to a specific state of the materials, hence it naturally applies to 
study mechanisms that are out of the scope of conventional poromechanics, such 
as phase transitions [8 or internal erosion 31 in saturated granular materials. The 
current formulation is restricted to quasi-static regimes. The extension to inertial 
regimes is currently in development. 
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